clear all;

# declare variables
filename = "e3registro";
report = "~/Documents/ITBA/ss/proyecto_ss/plots/report_E3";
clases = 7;

for i=1:clases
	cant(i) = 0;
	O(i) = 0;
	M(i) = 0;
	X(i) = 0;
endfor

#abro archivo para el report
report = fopen(report, "w");
fprintf(report, "Resultados para E3\n\n");

# abro, leo el archivo y lo cierro
file = fopen(filename);
fscanf(file, "%c", 66);	#consumo el header
i = 1;
time_max = 0;
time_min = 1000;
while (!feof(file))
	times(i) = fscanf(file, "%f", 1);
	if (times(i) > time_max)
		time_max = times(i);
	endif
	if (times(i) < time_min)
		time_min = times(i);
	endif
	i++;
endwhile
fclose(file);
n = i-1;


# computo las clases
dif = (time_max - time_min);
fprintf(report, "N\t\t%d\n", n);
fprintf(report, "Clases\t\t%d\n", clases);
fprintf(report, "Maximo\t\t%g\nMinimo\t\t%g\n", time_max, time_min);
for i=1:n
	for j=1:clases
		if (times(i) > ((j-1)/clases)*dif+time_min && times(i) < (j/clases)*dif+time_min)
			O(j) = O(j) + 1;
		endif
		M(j) = ((j-1)+0.5)*((time_max-time_min)/clases) + time_min;
	endfor
endfor


# computo los estadisticos de maxima verosimilitud
mu = 0;
sigma = 0;
for j=1:clases
	X(j) = M(j) * O(j);
	mu = mu + X(j);
endfor
mu = mu / n;
for j=1:clases
	sigma = sigma + (X(j) - mu)^2;
endfor
sigma = sqrt(sigma/n);

sigma = 0;
for i=1:n
sigma = sigma + ((times(i) - mu)^2);
endfor
sigma = sqrt(sigma/n);
fprintf(report, "mu\t\t%g\nsigma\t%g\n", mu, sigma);


# calculo la cantidad esperada para la normal
printf("j\tO(j)\tM(j)\t\tX(j)\t\ta\t\tb\t\tE(j)\n");
for j=1:clases
X(j) = O(j)*M(j);
a = ((j-1)/clases)*dif+time_min;
b = (j/clases)*dif+time_min;
E(j) = (normcdf(b,mu,sigma) - normcdf(a,mu,sigma)) * n;
printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\n",j,O(j),M(j),X(j),a,b,E(j));
endfor




# TEST CHI CUADRADO
# H0: PROVIENE DE UNA EXPONENCIAL
SumChi = 0;
fprintf(report, "\nTEST CHI CUADRADO\n");
fprintf(report, "H0: El tiempo entre arribos de E3 proviene de una distribucion normal N ~ (mu, sigma)\n");
fprintf(report, "Clases\tO(i)\tE(i)\t\tO(i)-E(i)\t(O(i)-E(i))^2\t\t(O(i)-E(i))^2 / E(i)\n");
for j=1:clases
	DifCuad(j) = ((O(j)-E(j))^2);
	Chi(j) = DifCuad(j) / E(j);
	SumChi = SumChi + Chi(j);
fprintf(report, "%d\t%d\t%g\t\t%g\t\t%g\t\t\t%g\n",j,O(j),E(j),O(j)-E(j),DifCuad(j),Chi(j));
endfor
fprintf(report, "\nEl estadistico Chi Cuadrado es %g\n",SumChi);
fprintf(report, "Buscar en tabla el Chi Critico para %d grados de libertad y el nivel de confianza deseado\n", clases-1);
						  
# TEST Q-Q
y = sort(times);
for i=1:n
	gamma(i) = (i-0.5)/n;
	gammaInv(i) = norminv(gamma(i), mu, sigma);
endfor
plot(gammaInv,y, '.');
legend("Plot Q-Q");
xlabel("Quantil Teorico");
ylabel("Quantil Observado");
print -deps "~/Documents/ITBA/ss/proyecto_ss/plots/PlotQQ_E3.eps";

# histograma
hist(times, clases);
legend("Histograma de tiempos entre llegadas al registro");
xlabel("Tiempo entre arribos");
ylabel("Frecuencia");
print -deps "~/Documents/ITBA/ss/proyecto_ss/plots/Histograma_E3.eps";

fclose(report);